#################################################################################### 
####################################################################################
####################################################################################
# Analyze HH-Level Data
####################################################################################
####################################################################################
library(here)
library(kableExtra)
library(lmtest)
library(multiwayvcov)
####################################################################################
####################################################################################
## Read in survey data (from: analysis_results.R, analysis_controls2.R)
dat = readRDS(here("data/dat_unrestricted.rds"))

####################################################################################
####################################################################################
## Table C.2 (two parts; requires post-processing)
tab = dat

## CHP Awareness Tables
cols = c("charity_for_profit","active_brac_lg")
tab[cols][is.na(tab[cols])] = -99

char_t = addmargins(table(tab[tab$treatment == 1,]$charity_for_profit), 1)
char_c = addmargins(table(tab[tab$treatment == 0,]$charity_for_profit), 1)

acti_t = addmargins(table(tab[tab$treatment == 1,]$active_brac_lg), 1)
acti_c = addmargins(table(tab[tab$treatment == 0,]$active_brac_lg), 1)

a_t = apply( round(prop.table(char_t[1:3]), digits = 2)*100, 1, function(u) paste0("(", sprintf( "%.0f%%", u ), ")") )
a_c = apply( round(prop.table(char_c[1:3]), digits = 2)*100, 1, function(u) paste0("(", sprintf( "%.0f%%", u ), ")") )
b_t = apply( round(prop.table(acti_t[1:3]), digits = 2)*100, 1, function(u) paste0("(", sprintf( "%.0f%%", u ), ")") )
b_c = apply( round(prop.table(acti_c[1:3]), digits = 2)*100, 1, function(u) paste0("(", sprintf( "%.0f%%", u ), ")") )

a_t[4] = ""
a_c[4] = ""
b_t[4] = ""
b_c[4] = ""

char = rbind(char_t,a_t, char_c, a_c)
acti = rbind(acti_t,b_t, acti_c, b_c)

colnames(char) = c("Don't Know", "Charity", "For-Profit", "Sum")
rownames(char) = c("Treatment", "", "Control", "")
colnames(acti) = c("Don't Know", "No", "Yes", "Sum")
rownames(acti) = c("Treatment", "", "Control", "")


stargazer(acti, type = "latex", label = "tab:mani_check_acti",
          #out = "tables/mani_check_acti.tex", 
          summary = F, title = "Respondents reporting active CHP in village")

## CHP Contact Tables
# 4 "Within the past one month" 3 "Within the past three months" 2 "Within the past six months" 1 "Within the past year" -1 "More than 1 year ago" 0 "Never" -99 "Don't Know"
ptvtme = median(tab[tab$treatment == 1,]$promoter_total_visits, na.rm = T)
ptvtmn = mean(tab[tab$treatment == 1,]$promoter_total_visits, na.rm = T)

ptvcme = median(tab[tab$treatment == 0,]$promoter_total_visits, na.rm = T)
ptvcmn = mean(tab[tab$treatment == 0,]$promoter_total_visits, na.rm = T)

ptv = cbind( rbind(ptvtme,ptvcme), rbind(ptvtmn,ptvcmn) )
rownames(ptv) = c("Treatment","Control")
colnames(ptv) = c("Median","Mean")

stargazer(ptv, type = "latex", label = "tab:mani_check_ptv",
          #out = "tables/mani_check_ptv.tex", 
          summary = F, title = "Care Received from CHP")

rm(acti, char, tab, ptv, ptvcme, ptvcmn, ptvtme, ptvtmn)

## CHP Satisfaction Plot (Treatment villages only)
chps <- dat[dat$treatment == 1 & !is.na(dat$promoter_satisfied),] %>% 
  group_by(promoter_satisfied) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count))

chp = ggplot(chps, aes(x=factor(promoter_satisfied),  y = perc*100)) + 
  geom_bar(stat="identity") +
  geom_vline(xintercept = mean(dat[dat$treatment == 1,]$promoter_satisfied, na.rm = T),
             color="blue", linetype="dashed", size=1) + 
  geom_text(data=chps, aes(label=paste0(round(perc*100,1),"%"),
                           y=perc*100-3), size=4, family="CM Roman") +
  geom_text(aes(x = (mean(dat[dat$treatment == 1,]$promoter_satisfied, na.rm = T)-.5), 
                label= paste0("mu =","=",round(mean(dat[dat$treatment == 1,]$promoter_satisfied, na.rm = T), 2)), y = 30), 
            colour="blue", angle=0, vjust = 0, parse = T,  size = 4, family="CM Roman") + 
  theme_bw() + 
  labs(title = NULL,
       x = NULL,
       y = NULL) + 
  scale_x_discrete(labels= c("(1)\nVery\ndissatisfied", "(2)\nSomewhat\ndissatisfied", "(3)\nNeither", 
                             "(4)\nSomewhat\n satisfied", "(5)\nVery\nsatisfied")) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(chps$perc*100)+1.5), labels=function(x) paste0(x,"%")) +
  theme(axis.title = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        plot.title = element_text(hjust = 0.5),
        text=element_text(family="CM Roman")) 

#ggsave(file = here("plots/chp_satisfaction",".pdf", sep = ""), chp, width = 5, height = 3)

rm(chp,chps)

####################################################################################
####################################################################################
## Manipulation  Checks

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Bushenyi","Ibanda","Jinja","Mbale","Mpigi","Pallisa","Mbarara","Shema",
         "village_distance_HF","village_distance_HOSP","village_distance_road","village_distance_electricity", "share2006")
dat[, paste0("c_", cols)] = lapply(dat[, cols], function(x) (x - mean(x))  )
covariates = c("treatment*c_village_distance_HF + treatment*c_village_distance_HOSP + treatment*c_village_distance_road + treatment*c_village_distance_electricity + treatment*c_sex + c_age*treatment + c_education_level*treatment + treatment*c_share2006 + c_Arua*treatment + 
               c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + 
               c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment","c_Arua*treatment + 
               c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + 
               c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment")


####################################################################################
####################################################################################
## Table C.3

## NGO Count (Total, Health, Index)
# Thinking about this list of NGOs that have been active in your village, did you or someone in your family receive services from any of these NGOs?
# 0 "No" 1 "Once or twice" 2 "More than twice" 3 "More than five times" 4 "More than ten times"
for (i in names(dat[,c("ngo_count_health", "ngo_count_nonhealth", "ngo_contact_health", "ngo_contact_nonhealth")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

stargazer(m.ngo_count_health.c, m.ngo_count_health.t, m.ngo_contact_health.c, m.ngo_contact_health.t,  
          m.ngo_count_nonhealth.c, m.ngo_count_nonhealth.t, m.ngo_contact_nonhealth.c, m.ngo_contact_nonhealth.t, 
          
          se = list(c.ngo_count_health.c[,2], c.ngo_count_health.t[,2], c.ngo_contact_health.c[,2], c.ngo_contact_health.t[,2], 
                    c.ngo_count_nonhealth.c[,2], c.ngo_count_nonhealth.t[,2], c.ngo_contact_nonhealth.c[,2], c.ngo_contact_nonhealth.t[,2] ),
          
          p = list(c.ngo_count_health.c[,4], c.ngo_count_health.t[,4], c.ngo_contact_health.c[,4], c.ngo_contact_health.t[,4], 
                   c.ngo_count_nonhealth.c[,4], c.ngo_count_nonhealth.t[,4], c.ngo_contact_nonhealth.c[,4], c.ngo_contact_nonhealth.t[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", here("tables/manipulation_1.tex"),
          label = "tab:manipulation_1",column.sep.width = "1pt",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of CHP Intervention on Perceptions of NGO Activity",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("NGO Count","NGO Contact","NGO Count","NGO Contact"), column.separate = c(2,2,2,2),
          add.lines = list(c("Covariates", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")))

# Post-processing: manipulation_1.tex
# & \multicolumn{4}{c}{Health NGOs} & \multicolumn{4}{c}{Non-Health NGOs} \\ 
# \cmidrule(lr){2-5} \cmidrule(l){6-9} 

rm(list =  grep("^(c|m)\\.", ls(), value = T))

####################################################################################
####################################################################################
## Table C.4

## Personal and Community Benefits (benefits from NGOs... Very big - Not big at all)
for (i in names(dat[,c("personal_benefits","community_benefits")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

stargazer(m.personal_benefits.c, m.personal_benefits.t, m.community_benefits.c, m.community_benefits.t,  
          
          se = list(c.personal_benefits.c[,2], c.personal_benefits.t[,2], c.community_benefits.c[,2], c.community_benefits.t[,2] ),
          
          p = list(c.personal_benefits.c[,4], c.personal_benefits.t[,4], c.community_benefits.c[,4], c.community_benefits.t[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = here("tables/manipulation_2.tex"),
          label = "tab:manipulation_2", column.sep.width = "1pt",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of CHP Intervention on Perceptions of Benefits from NGOs",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("Respondent's Household","Respondent's Community"), column.separate = c(2,2),
          add.lines = list(c("Covariates", "No", "Yes", "No", "Yes")))

rm(list =  grep("^(c|m)\\.", ls(), value = T))
####################################################################################